{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interoperability with Atomic Simulation Environment\n",
    "\n",
    "\\[[View with IPython nbviewer](http://nbviewer.ipython.org/github/libAtoms/QUIP/blob/public/doc/Tutorials/quippy-ase-interoperability.ipynb)\\]\n",
    "\n",
    "- quippy now uses the standard `ase.atoms.Atoms` class to represent Atomic structures\n",
    "- quippy `Potential` objects can be used as ASE calculators, and vice-versa\n",
    "- Can use standard ASE tools, plus communicate with other packages using ASE as *lingua franca*\n",
    "\n",
    "## Example: vacancy formation energy\n",
    "\n",
    "- Generate structure with `ASE` lattice tools\n",
    "- Stillinger-Weber potential implementation from `QUIP`\n",
    "- Elastic constant fitting routine from `matscipy`, internal relaxations with `ASE` FIRE minimiser"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n",
      "      Step     Time          Energy         fmax\n",
      "BFGS:    0 21:22:10      -34.635777        0.3247\n",
      "BFGS:    1 21:22:10      -34.644590        0.1504\n",
      "BFGS:    2 21:22:10      -34.646997        0.0001\n",
      "      Step     Time          Energy         fmax\n",
      "BFGS:    0 21:22:10      -34.670667        0.1584\n",
      "BFGS:    1 21:22:10      -34.672777        0.0749\n",
      "BFGS:    2 21:22:10      -34.673385        0.0000\n",
      "      Step     Time          Energy         fmax\n",
      "BFGS:    0 21:22:10      -34.678737        0.0000\n",
      "      Step     Time          Energy         fmax\n",
      "BFGS:    0 21:22:10      -34.660845        0.1508\n",
      "BFGS:    1 21:22:10      -34.662784        0.0742\n",
      "BFGS:    2 21:22:10      -34.663403        0.0000\n",
      "      Step     Time          Energy         fmax\n",
      "BFGS:    0 21:22:10      -34.617822        0.2945\n",
      "BFGS:    1 21:22:10      -34.625262        0.1477\n",
      "BFGS:    2 21:22:10      -34.627763        0.0001\n",
      "Fitting C_11\n",
      "Strain array([-0.02, -0.01,  0.  ,  0.01,  0.02])\n",
      "Stress array([-2.56044687, -1.01247671,  0.5027424 ,  1.98366491,  3.42893711]) GPa\n",
      "Cij (gradient) / GPa    :     149.74909575669915\n",
      "Error in Cij / GPa      :     1.1696996170085603\n",
      "Correlation coefficient :     0.9999084935045536\n",
      "Setting C11 (1) to 0.934660 +/- 0.007301\n",
      "\n",
      "\n",
      "Fitting C_21\n",
      "Strain array([-0.02, -0.01,  0.  ,  0.01,  0.02])\n",
      "Stress array([-1.07663577, -0.26643655,  0.5027424 ,  1.23345414,  1.92818198]) GPa\n",
      "Cij (gradient) / GPa    :     75.0952617697293\n",
      "Error in Cij / GPa      :     1.3149075235415504\n",
      "Correlation coefficient :     0.9995404242109733\n",
      "Setting C21 (7) to 0.468708 +/- 0.008207\n",
      "\n",
      "\n",
      "Fitting C_31\n",
      "Strain array([-0.02, -0.01,  0.  ,  0.01,  0.02])\n",
      "Stress array([-1.07663577, -0.26643655,  0.5027424 ,  1.23345414,  1.92818198]) GPa\n",
      "Cij (gradient) / GPa    :     75.09526176972929\n",
      "Error in Cij / GPa      :     1.31490752354155\n",
      "Correlation coefficient :     0.9995404242109733\n",
      "Updating C31 (7) with value 0.468708 +/- 0.008207\n",
      "\n",
      "\n",
      "Fitting C_44\n",
      "Strain array([-0.02, -0.01,  0.  ,  0.01,  0.02])\n",
      "Stress array([-1.13572340e+00, -5.65842409e-01, -9.46072689e-15,  5.60142655e-01,\n",
      "        1.11304586e+00]) GPa\n",
      "Cij (gradient) / GPa    :     56.23523568430933\n",
      "Error in Cij / GPa      :     0.19437884854805132\n",
      "Correlation coefficient :     0.9999820790695022\n",
      "Setting C44 (4) to 0.350993 +/- 0.001213\n",
      "\n",
      "\n",
      "[[b C11 b C12 b C12 b     b     b    ]\n",
      " [b C12 b C11 b C12 b     b     b    ]\n",
      " [b C12 b C12 b C11 b     b     b    ]\n",
      " [b     b     b     b C44 b     b    ]\n",
      " [b     b     b     b     b C44 b    ]\n",
      " [b     b     b     b     b     b C44]]\n",
      "\n",
      " = \n",
      "\n",
      "[[149.75  75.1   75.1    0.     0.     0.  ]\n",
      " [ 75.1  149.75  75.1    0.     0.     0.  ]\n",
      " [ 75.1   75.1  149.75   0.     0.     0.  ]\n",
      " [  0.     0.     0.    56.24   0.     0.  ]\n",
      " [  0.     0.     0.     0.    56.24   0.  ]\n",
      " [  0.     0.     0.     0.     0.    56.24]]\n",
      "C_11 = 149.75 +/- 1.17 GPa\n",
      "C_12 = 75.10 +/- 1.31 GPa\n",
      "C_44 = 56.24 +/- 0.19 GPa\n",
      "PreconLBFGS:   0  21:22:10     -927.087471       0.8332\n",
      "estimate_mu(): mu=2.3315998829549316, mu_c=1.0\n",
      "PreconLBFGS:   1  21:22:10     -927.611735       0.1487\n",
      "PreconLBFGS:   2  21:22:11     -927.643943       0.0642\n",
      "PreconLBFGS:   3  21:22:11     -927.652380       0.0393\n",
      "PreconLBFGS:   4  21:22:11     -927.655827       0.0210\n",
      "PreconLBFGS:   5  21:22:11     -927.656988       0.0092\n",
      "PreconLBFGS:   6  21:22:11     -927.657131       0.0052\n",
      "PreconLBFGS:   7  21:22:11     -927.657193       0.0030\n",
      "PreconLBFGS:   8  21:22:11     -927.657220       0.0011\n",
      "PreconLBFGS:   9  21:22:11     -927.657224       0.0009\n",
      "PreconLBFGS:  10  21:22:11     -927.657226       0.0004\n",
      "PreconLBFGS:  11  21:22:11     -927.657226       0.0002\n",
      "PreconLBFGS:  12  21:22:11     -927.657226       0.0001\n",
      "PreconLBFGS:  13  21:22:11     -927.657226       0.0001\n",
      "PreconLBFGS:  14  21:22:11     -927.657226       0.0000\n",
      "PreconLBFGS:  15  21:22:11     -927.657226       0.0000\n",
      "PreconLBFGS:  16  21:22:11     -927.657226       0.0000\n",
      "PreconLBFGS:  17  21:22:11     -927.657226       0.0000\n",
      "PreconLBFGS:  18  21:22:12     -927.657226       0.0000\n",
      "PreconLBFGS:  19  21:22:12     -927.657226       0.0000\n",
      "PreconLBFGS:  20  21:22:12     -927.657226       0.0000\n",
      "PreconLBFGS:  21  21:22:12     -927.657226       0.0000\n",
      "SW vacancy formation energy 4.33384020178346 eV\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "from ase.build import bulk\n",
    "from ase.optimize import BFGS\n",
    "from ase.optimize.precon import PreconLBFGS\n",
    "from quippy.potential import Potential\n",
    "\n",
    "si = bulk('Si', a=5.44, cubic=True)\n",
    "sw_pot = Potential('IP SW', \n",
    "                   param_filename='../../share/Parameters/ip.parms.SW.xml')\n",
    "si.set_calculator(sw_pot)\n",
    "e_bulk_per_atom = si.get_potential_energy()/len(si)\n",
    "\n",
    "# call general purpose elastic constants calculator \n",
    "#   using ASE Atoms and QUIP Potential\n",
    "from matscipy.elasticity import fit_elastic_constants\n",
    "Cij = fit_elastic_constants(si, optimizer=BFGS,\n",
    "                            symmetry='cubic', logfile='-')\n",
    "vac1 = si.copy()\n",
    "vac1 *= (3, 3, 3)\n",
    "half_cell = np.diag(vac1.cell)/2.\n",
    "vac_atom = ((vac1.positions - half_cell)**2).sum(axis=1).argmin()\n",
    "del vac1[vac_atom]\n",
    "\n",
    "vac1.set_calculator(sw_pot)\n",
    "vac1.rattle(0.01)\n",
    "opt = PreconLBFGS(vac1) # big cell, use preconditioned minimiser\n",
    "opt.run(fmax=1e-6)\n",
    "e_vac = vac1.get_potential_energy() - e_bulk_per_atom*len(vac1)\n",
    "print('SW vacancy formation energy', e_vac, 'eV')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
